function rt_2dhist(axesH, x, y, xbins, ybins, levels)
    
    [rows, cols] = size(x);
    if cols ~= 1
        x = x';
    end
    [rows, cols] = size(y);
    if cols ~= 1
        y = y';
    end
    n          = length(y);
    mYX        = [y, x];

    vYEdge     = linspace(ybins(1), ybins(2), ybins(3));
    vXEdge     = linspace(xbins(1), xbins(2), xbins(3));
    mHist2d    = hist2d(mYX, vYEdge, vXEdge);
    maxHist    = max(max(mHist2d));
    contourint = 1:1:maxHist;

    nXBins     = length(vXEdge);
    nYBins     = length(vYEdge);
    vXLabel    = 0.5*(vXEdge(1:(nXBins-1)) + vXEdge(2:nXBins));
    vYLabel    = 0.5*(vYEdge(1:(nYBins-1)) + vYEdge(2:nYBins));

    mHist2d_prob = zeros(size(mHist2d));
    for i = 1:length(contourint)
        tmp = find(mHist2d >= contourint(i));
        mHist2d_prob(tmp)= 1 - (sum(mHist2d(tmp))/n);
    end

    ind_tmp      = find(mHist2d_prob >= 0);
    [X, Y]        = meshgrid(vXLabel, vYLabel);
    [Xsm, Ysm]    = meshgrid(linspace(vXLabel(1), vXLabel(end), 1024), ...
                             linspace(vYLabel(1), vYLabel(end), 1024));
    mHist2d_prob = griddata(X(ind_tmp), Y(ind_tmp),...
                            mHist2d_prob(ind_tmp), Xsm, Ysm, 'cubic');
    [~, h] = contour(axesH, Xsm, Ysm, mHist2d_prob, 0.05:0.05:0.95, 'fill', 'on');
    CM = get(h,'ContourMatrix');
    CM = CM(:,(CM(1,:) > min(x)) & (CM(1,:) < max(x)));
    xdat = [min(CM(1,:)) max(CM(1,:))];
    ydat = [min(CM(2,:)) max(CM(2,:))];
    set(gca(), 'XLim', fix_lims(get(gca, 'XLim'), xdat));
    set(gca(), 'YLim', fix_lims(get(gca, 'YLim'), ydat));
    caxis([0 1]);
    colormap(jet);
    contour(axesH, Xsm, Ysm, mHist2d_prob, levels, 'color', 'k');
end

function newlims = fix_lims(oldlims, data)
    newlims = oldlims;
    range = oldlims(2) - oldlims(1);
    % fix lower lim, if necessary
    if data(1) - oldlims(1) > 0.1*range
        newlims(1) = data(1) - 0.1*range;
    end
    if oldlims(2) - data(2) > 0.1*range
        newlims(2) = data(2) + 0.1*range;
    end
end